subroutine boundary_check(px,py,pz,removed)
  use param
  use box_param,only:bwidth
  use box_param,only:bound
  implicit none

  real(dbl),intent(inout) :: px,py,pz
  logical,intent(out) :: removed

  removed=.false.

  if (px.gt.0.and.px.le.bwidth(1)) then
  else if (px.le.0) then
     ! left in x
     if (bound(1).eq.0) then
        ! Free boundary
        removed=.true.
        return
     else
        ! Periodic boundary
        do
           px=px+bwidth(1)
           if (px.gt.0) exit
        end do
     end if
  else
     ! right in x
     if (bound(1).eq.0) then
        ! Free boundary
        removed=.true.
        return
     else
        ! Periodic boundary
        do
           px=px-bwidth(1)
           if (px.gt.0.and.px.le.bwidth(1)) exit
        end do
     end if
  end if

  if (py.gt.0.and.py.le.bwidth(2)) then
  else if (py.le.0) then
     ! left in y
     if (bound(2).eq.0) then
        ! Free boundary
        removed=.true.
        return
     else
        ! Periodic boundary
        do
           py=py+bwidth(2)
           if (py.gt.0) exit
        end do
     end if
  else
     ! right in y
     if (bound(2).eq.0) then
        ! Free boundary
        removed=.true.
        return
     else
        ! Periodic boundary
        do
           py=py-bwidth(2)
           if (py.gt.0.and.py.le.bwidth(2)) exit
        end do
     end if
  end if

  if (pz.gt.0.and.pz.le.bwidth(3)) then
  else if (pz.le.0) then
     ! left in z
     if (bound(3).eq.0) then
        ! Free boundary
        removed=.true.
        return
     else
        ! Periodic boundary
        do
           pz=pz+bwidth(3)
           if (pz.gt.0) exit
        end do
     end if
  else
     ! right in z
     if (bound(3).eq.0) then
        ! Free boundary
        removed=.true.
        return
     else
        ! Periodic boundary
        do
           pz=pz-bwidth(3)
           if (pz.gt.0.and.pz.le.bwidth(3)) exit
        end do
     end if
  end if

  return

end subroutine boundary_check
